Treating instabilities in a hyperbolic formulation of Einstein's equations 
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We have recently constructed a numerical code that evolves a spherically symmetric spacetime 
using a hyperbolic formulation of Einstein's equations. For the case of a Schwarzschild black hole, 
this code works well at early times, but quickly becomes inaccurate on a time scale of 10 — lOOAf, 
where M is the mass of the hole. We present an analytic method that facilitates the detection 
of instabilities. Using this method, we identify a term in the evolution equations that leads to a 
rapidly-growing mode in the solution. After eliminating this term from the evolution equations 
by means of algebraic constraints, we can achieve free evolution for times exceeding 10 000M. We 
discuss the implications for three-dimensional simulations. 
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I. INTRODUCTION 

When solving Einstein's equations as an initial value 
problem, one considers spacetime as a foliation of space- 
like hypersurfaces, or "slices" . Einstein's equations then 
separate into two types: constraint equations, which re- 
late the dynamical variables on each particular slice, and 
evolution equations, which describe how these variables 
propagate from one slice to the next. The constraints 
are analogous to the divergence equations in Maxwell's 
theory, and the evolution equations are analogous to the 
curl equations. 

As in Maxwell's theory, the evolution equations admit 
solutions that violate the constraints. However, if the 
constraints are satisfied on the initial slice and on all spa- 
tial boundaries, then the evolution equations guarantee 
that the constraints are satisfied elsewhere. This permits 
numerical solution schemes in which only the evolution 
equations are explicitly solved at each time step. 

Such "free evolution" schemes are desirable for several 
reasons. First, the constraints are typically nonlinear el- 
liptic equations, which are difficult and costly to solve on 
a computer, especially in the general case of three spatial 
dimensions. Second, a free evolution scheme allows one 
to track numerical errors by monitoring the constraints 
at each time step. 

For numerical evolution of black holes, an additional 
advantage of a free evolution scheme is that one can, 
in principle, excise a black hole from the spacetime and 
evolve only the exterior region, and one can do so with- 
out imposing explicit boundary conditions on the hori- 
zon. This is the basis for so-called "apparent horizon 
boundary condition" methods, which are thought to be 
crucial for long-term numerical evolution of black hole 
spacetimes [DHI. However, excising a black hole from 
a spacetime is known to be mathematically well-defined 



only if the evolution equations are hyperbolic and if the 
characteristic curves of the hyperbolic system are "phys- 
ical", that is, if they lie within the local light cone. In 
this case, the structure of the equations guarantees that 
no information, including gauge information, can emerge 
from the hole. For non-hyperbolic representations of gen- 
eral relativity such as the usual ADM |HJ formulation, 
the evolution equations are of no mathematical type for 
which well-posedness has been proven, so the suitabil- 
ity of these formulations for black hole excision must be 
determined empirically on a case-by-case basis. It is in 
part for this reason that much attention has been recently 
focused on hyperbolic representations of Einstein's equa- 
tions 

A key stumbling block in numerical work, particularly 
in finite-difference solutions of initial value problems, is 
the tendency for numerical computations to become un- 
stable. Instabilities have many origins, and the cause 
of any particular instability found in a numerical code 
is often difficult to deduce. Furthermore, if the desired 
analytic solution is unknown, it can be difficult to dis- 
tinguish between an instability and a case in which the 
analytic solution simply grows without bound. Exam- 
ples of the latter include systems that evolve to physical 
singularities (e.g., Oppenheimer-Snyder collapse evolved 
using geodesic slicing) and those that evolve toward co- 
ordinate singularities (e.g., a Schwarzschild black hole 
evolved with maximal time slicing, and several harmonic- 
slicing examples that become singular for certain choices 
of the initial lapse function |2(],|2l| ) . When diagnosing in- 
stabilities in numerical simulations, it is therefore prefer- 
able to study instances in which the analytic solution is 
known and well-behaved. 

We distinguish between two types of instabilities: a 
type in which the numerical finite-difference equations 
admit rapidly-growing solutions that do not satisfy the 
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underlying continuum differential equations, and a type 
in which the continuum equations themselves admit 
growing modes that are absent in the desired solution 
but are excited by numerical perturbations. An example 
of the former type, which we will call a numerical instabil- 
ity because it depends on the numerical finite-difference 
equations, is the well-known Courant instability that can 
arise in explicit finite-difference solutions of hyperbolic 
PDEs. The high-frequency modes that characterize a 
Courant instability do not satisfy the underlying differ- 
ential equations. 

The latter type, which we will call a "continuum" in- 
stability because the unstable mode satisfies the contin- 
uum differential equations, commonly occurs in systems 
of equations that admit both well-behaved and grow- 
ing solutions. Although one might be interested in the 
well-behaved solution, the growing mode eventually dom- 
inates if it at any time acquires a nonzero amplitude 
via numerical errors. A simple example is the equa- 
tion y — y/9 with initial conditions y = 1 , y = — 1/3. 
For these initial conditions the unique analytic solution 
is y = e~*/ 3 , but a naive numerical integration of this 
problem is unstable as it proceeds forward in time be- 
cause numerical perturbations excite the growing solu- 
tion y = e*/ 3 . 

For numerical solutions of Einstein's equations, a con- 
tinuum instability may be due to a gauge mode excited by 
inaccuracies in numerically-determined coordinate condi- 
tions. Or, in the case of a free evolution scheme, it may 
be caused by a rapidly-growing mode that satisfies the 
evolution equations but violates the constraints. This 
latter case is possible despite the fact that the evolu- 
tion equations preserve the constraints, because in nu- 
merical computations neither the evolution equations nor 
the constraints are exactly satisfied. Constraint-violating 
modes have been discussed in the literature p2|-p4| but 
their importance for numerical free evolution schemes re- 
mains controversial. 

Eliminating a continuum instability often requires a 
different approach than removing a numerical one, be- 
cause these two types of instability stem from quite dif- 
ferent sources. To remove a numerical instability, one 
must change the numerical algorithm (or details of the 
algorithm such as the size of the time step) that is used 
to solve the equations, so that this algorithm no longer 
introduces growing modes. To remove a continuum in- 
stability, one must either remove the numerical pertur- 
bations that excite the undesired solution of the contin- 
uum equations, change the numerical scheme in order to 
damp out this solution, or modify the continuum equa- 
tions themselves (possibly including the choice of gauge) 
so that no growing solution is present. 

In this paper we examine instabilities in a numerical 
free evolution code that solves a spherically symmetric 
black-hole spacetime. Our code, which has been de- 
scribed in detail elsewhere 0, is based on a hyperbolic 
formulation of general relativity (the "Einstein-Ricci" 
or "ER" formulation) originally proposed by Choquet- 



Bruhat and York ,[L3|] . For short integration times our 
code performs well, but we show in section [II that for 
the case of a Schwarzschild black hole it becomes unsta- 
ble and terminates on a time scale of 10 — 100M, where 
M is the mass of the hole. This occurs even in a gauge in 
which the analytic solution is regular at the horizon and 
time- independent. The rate at which our errors grow is 
independent of the numerical time discretization At and 
the spatial discretization Ar, suggesting that the growth 
is due to a continuum instability rather than a numerical 
one. 



In section IV we present a method of analyzing the 
evolution equations that facilitates the detection of con- 
tinuum instabilities. In the simplest application of this 
method we consider each ER evolution equation sepa- 
rately. For each equation, we examine the free evolu- 
tion of the ER variable governed by that equation, treat- 
ing all other ER variables as fixed and given by the 
Schwarzschild solution. We ask whether perturbations 
of the evolved ER variable about its Schwarzschild value 
grow rapidly with time. We find that most of the ER 
equations, when treated individually in this manner, are 
stable, but that one of the ER equations is sensitive to 
a continuum instability. A single term on the right-hand 
side of the unstable equation is responsible for the grow- 
ing mode. 

In section ^ we construct a modified set of evolu- 
tion equations that no longer contain this troublesome 
term. This is done primarily by using algebraic con- 
straints to rewrite the right-hand side of one equation. 
We find that numerical free evolution of the modified 
set of equations remains accurate for times in excess of 
lOOOOAf. This substantial improvement indicates that 
the rapidly-growing mode found by our analysis in sec- 
tion IV is the dominant instability afflicting free evolu- 
tion of the unmodified ER equations. In section VI we 



discuss our method of stability analysis and apply it to 
the three-dimensional Einstein-Ricci equations, as well as 
to the Einstein-Bianchi [[l9| and ADM systems. We dis- 
cuss the implications for three-dimensional free evolution 
schemes. 



II. EQUATIONS 
A. The ER Formalism 

Here we summarize the fundamental variables and 
equations used in the ER representation of general rela- 
tivity. For details of the ER formulation and a derivation 
of the equations, see [ ^2| , p^[ . 

We write the metric in the usual 3+1 form 



ds 2 



N 2 dt A + giJ {dx l + f3 % dt){dx 3 + P>dt), (1) 



where N is the lapse function, (3 l is the shift vector, and 
gij is the three-metric on a spatial hypersurface of con- 
stant t. 
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Define the variables 





= —-N 1 dnqi4, 


(2a) 






(2b) 


M kij 


= D k K i: j, 


(2c) 


a,i 


= A(lniV), 


(2d) 




= N^doGi, 


(2e) 


a jj 




(2f) 



Here D is the three-dimensional covariant derivative com- 
patible with the three-metric gij , the time derivative op- 
erator is 



d = 



d_ 

dt 



£ 



/3. 



(3) 



and £ denotes a Lie derivative. The quantity is the 
usual extrinsic curvature. 

The vacuum evolution equations for the general three- 
dimensional case can be found in |l2|Jl3| j7[|. The vacuum 
constraint equations include 



Rij ^ij 



ah 



3i 

aoi + Ha 



M - J 

1V1 IJ ; 



(4a) 
(4b) 
(4c) 
(4d) 



where Rij is the three-dimensional Ricci tensor formed 



from the three-metric g^. Eqs. (4a)-(|4c[) follow from 
the Gauss-Codazzi-Ricci equations for embeddinga foli- 
ation into a higher-dimensional space, and Eq. (Eq) fol- 
lows from harmonic time slicing. Additional constraints 
that must be satisfied at all times are the definitions (pq 1 ), 
(pc|), and (^|), and the usual relation between T k ij and 
derivatives of g^ . 



B. Spherical Symmetry 

The spherically symmetric three-metric can be written 
in the general form 



< 3 >ds 2 -- 


= A 2 dr 2 + B 2 r 2 {d9 2 + sin 2 6d4> 2 ), 


(5) 


(r,M) 


are the usual spherical coordinates. 


Define 


T r T = 


2BrT e 9r = 2BrT% r 






2 A 2 2 A 2 
Br Tee - Br sin 2 / **' 


(6a) 


ax 


6 <h 

a e = a <j>, 


(6b) 


I . \ = 


Li $ = -L 0, 


(6c) 


K T = 


K% = K*+, 


(6d) 


M rT = 


M r \ = M/ 0) 


(6e) 


M Tr = 


M% r = M+fr, 


(6f) 



where the subscript T denotes "transverse" . 

The evolution equations can be written in the form 



d A = -NAK r r , 
d Br = -NBrK T , 



d K rr 
d K T 
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do<l r 

d a T 
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T r T ( L rr T 
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d 
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d r CLQr 



(7a) 
(7b) 
(7c) 
(7d) 
(7e) 
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(7k) 
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IVr 

rr Br 
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\ Br 

d 

d Mrrr = N—L rr + N \(a r - 2Y r rr )L„ 
or 



(71) 



(7m) 



+ 2Kr{K rr a r + Mrrr)] , 

N d 

d Q L rr = — —Mrrr + N \L r r{^K T - §KZ) + 8M rT a 
A z or 



+ 



do M r T = N—L T 
or 



3(o r - T r rr ) + 
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(7n) 
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d * L T=-£Q- r MrT + N 



LtK 
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1 



-(Lrr + a rr ){K r r - K T )^ - 2K T 3 
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+2a T {K r r + K T ) - 2K r r 2 K T + Kf 
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rT 
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T t t + T r rr T rT 



2A 2 Br 



K rr {2K T -K r r )-L rr = Q, 

2" 



-_r rT + r rr r rT - — 



1 



B 2 r 2 
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2Lt + -j? + 2A> 2 + iC 2 + 2a T 
A z 

1 



0, 



0, 



M rT - M Tr = 0, 



A 2 



2M rT = 0. 



(8a) 

(8b) 

(8c) 
(8d) 
(8e) 



The additional constraints (|2cJ), (2d), (|2J), and the usual 
relation between T a and derivatives of q,,, take the form 
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9ij 
0. 
0, 
0, 
0, 
0, 
0. 
0, 
0. 



(8f) 

(8g) 
(8h) 
(8i) 

(8j) 
(8k) 
(81) 
(8m) 



III. FREE EVOLUTION OF ER SYSTEM 
A. Method 

We solve the spherically symmetric ER evolution equa- 
tions (|?]) at every time step using the causal differencing 
method described in . The constraints are satisfied on 
the initial time slice but are not solved explicitly during 
the evolution. 

The inner boundary of the numerical domain is a sur- 
face that remains within a grid spacing of the apparent 



horizon, r = rAH- Because the apparent horizon is an 
outgoing null or spacelike surface, the hyperbolic evolu- 
tion equations require no boundary condition there. The 
outer boundary is an arbitrary spherical surface far from 
the black hole at r = r max . At the outer boundary, 
we use the "extended Robin" condition discussed in [[7|. 
This outer boundary condition does not properly handle 
wavelike behavior, but in practice it is adequate for the 
cases shown here. 

The lapse function can be freely specified on the initial 
time slice, and is subsequently determined by the har- 
monic time slicing condition Ot = 0. The shift is chosen 
to satisfy the minimal strain equation p5[ | . This equation 
minimizes the average change in the three-metric as one 
evolves from one time slice to the next, and is used to 
provide a shift vector that does not produce coordinate 
singularities. The minimal strain equation requires two 
boundary conditions, for which we choose 



N 
~A 







at r = r A H, 



dr 



(r 2 (3 r ) =0 at r = r n 



(9) 
(10) 



The inner boundary condition ensures that at the appar- 
ent horizon, the coordinates move outward at the local 
speed of light, c = N/A. This prevents the coordinates 
from falling into the black hole. The outer boundary con- 
dition ensures that the shift falls off like r~ 2 , in accor- 
dance with the time-independent Schwarzschild solution 
written in harmonic slicing (Eqs. (11) below). We use a 
feedback technique M to keep the horizon near r = 2M. 



B. Initial data 

Our initial data are chosen on a time slice correspond- 
ing to a well-behaved, fully time-independent harmonic 
foliation of the Schwarzschild geometry (cf. refs. [p6|-p8|). 
Such a slice penetrates the event horizon without encoun- 
tering a coordinate singularity, and extends to the phys- 
ical singularity at r = 0. With an appropriate choice of 
spatial coordinates on the slice, all dynamical variables 
are time-independent [E8| and are given by 
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(111) 
(11m) 

(11*) 
(llo) 

(HP) 
(llq) 



where M is the mass of the black hole. One can ex- 
plicitly check the time-independence of this solution by 
inserting ( |ll| ) into the ER evolution equations (0) and 
verifying that all time derivatives are zero. Note that ( O ) 
satisfies the minimal strain shift condition, as does any 
time- independent solution of Einstein's equations. 



C. Results 



Figure ^ shows the error in the metric function A as a 
function of time. We plot the quantity \ A — A 811 ^, where 



A an is the analytic value of A given by (|lla| ), and the £2 
norm of a quantity q is defined by 



N 



(12) 



The sum is over all grid points that contain valid data 
(i.e., all grid points outside the horizon). The quantity 
\A — yl an |2 is shown for several different grid resolutions, 
each with the same Courant factor At/Ar. 

At early times, the error in A varies with resolution 
like O(Ar) 2 , as expected for our second-order convergent 




FIG. 1. The £2 norm of A minus its analytic solution (11a), 
shown as a function of time for five grid resolutions. The outer 
boundary is at r max = 64M and the Courant factor At/Ar 
is 3/4. All five plots terminate when the code crashes. 



numerical method. However, after about 10-30M the 
error grows rapidly, approximately like t at late times. 
The growth rate is independent of the grid resolution. 
Eventually, when errors have become sufficiently large, 
the code crashes, typically because it fails to locate an 
apparent horizon. 

It is common for numerical finite-difference schemes to 
produce solutions with errors that grow as truncation er- 
ror accumulates. However, such growth is typically linear 
in time, with a slope proportional to (At) 2 (for a second- 
order scheme), and can be easily defeated by increasing 
the resolution. In contrast, Figure [l] shows a more rapid 
growth rate that increases with time, indicating that we 
are observing something other than accumulating trun- 
cation errors. 

In Figure ^| we plot the error in A as a function of 
radius for several different times. The error is greatest 
near the horizon and remains smooth in both space and 
time as it grows. The fact that our errors are largest 
near the black hole does not necessarily indicate that 
the instability is somehow associated with our treatment 
of the inner boundary; one expects numerical errors to 
be greater for smaller values of r simply because most 
quantities in (11) behave like 1 jr n with positive n. 

Other quantities behave much like the error in A. In 
Figure || we plot the error in L rr with respect to the ana- 
lytic solution (111), and in Figure ^ we plot the left-hand 
side of the Hamiltonian constraint (|k]). Both quantities 
are approximately second-order convergent, but at late 
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at a rate independent of the grid resolution. 
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FIG. 2. Error in A as a function of coordinate radius, for 
the Ar/M = 1/16 case shown in Figure |l[ The function 
A — A an is plotted at five times. The error grows rapidly but 
smoothly until the code crashes. 
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FIG. 3. The £2 norm of L rr minus its analytic solution ( |llj ) 
as a function of time, shown for the same cases as plotted in 
Figure |l| 



times they increase rapidly (faster than linearly) in time 
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FIG. 4. The £2 norm of the Hamiltonian constraint (pq) 
versus time, shown for the same five cases as in Figure jjj. 

Figures suggest that the instability is not purely 
numerical. Numerical instabilities typically grow like e n , 
where n is the number of time steps. Consequently, for a 
numerical instability one expects that reducing the time 
discretization At would make the instability grow faster 
as a function of time, because integrating to a particular 
value of t requires more steps. However, in Figures |l[- 
||, At is decreased with each finer grid resolution, but 
the growth rate is unaffected. Similarly, at late times 
we see no change in the growth rate if we vary At while 
keeping the grid resolution fixed, as shown in Figure |^. 
Instead, for At — > our errors converge to a limit (this 
is simply the limit in which numerical truncation error is 
dominated by Ar instead of At). 

Our results instead suggest that our code suffers from 
a continuum instability. In this case, the code should re- 
main second-order convergent and the growth rate of er- 
rors should depend only on the continuum equations and 
not on numerical parameters like Ar or At. A smaller 
At or Ar should not intensify the instability, but instead 
should improve our simulations by virtue of reducing the 
numerical perturbations that excite the offending mode. 
Our results are consistent with these expectations. 

One possible source of a continuum instability is 
a rapidly-increasing constraint-violating solution of the 
evolution equations that is being excited by numerical 
perturbations. Another is a gauge mode that is not 
present in the analytic solution. In the case of a gauge 
mode, one would expect gauge-invariant quantities to re- 
main relatively unaffected while other quantities blow up. 



G 



0.01 - 



0.001 =- 



< 
< 




0.0001 : 



FIG. 5. The 1% norm of A — A an versus time shown for five 
different values of At, each with Ar/M = 1/16. The outer 
boundary is at r — 64M. At late times, the dominant error 
is independent of At. 



However, at late times, both gauge-dependent quantities 
like L rr (Figure |3|) and gauge-invariant quantities like the 
Hamiltonian constraint (Figure |J) increase rapidly with 
time at approximately the same rate. 



IV. STABILITY OF INDIVIDUAL EVOLUTION 
EQUATIONS 

To gain further insight into the nature of the instabil- 
ity, we consider each ER evolution equation separately. 
For each evolution equation, we treat the ER variable 
governed by that equation as freely evolving, but we fix 
the remaining ER variables to the analytic expressions 
given in ([ll]). In this way we can study the stability 
of each individual evolution equation in the absence of 
all couplings to other equations. Although this analysis 
will not shed light on any instabilities that are caused 
by these couplings, it is likely that if any of the evolution 
equations are found to be unstable individually, they will 
remain unstable when coupled to the other equations. 

We note that the method of analysis described below 
can also be used to examine coupled sets of equations as 
long as the couplings do not arise from derivative terms — 
this is described in more detail in Appendix |a|. However, 
we will see that treating one equation at a time is suffi- 
cient for the case discussed here. 

Let y represent any of the ER variables that evolve ac- 
cording to (0). If all ER variables other than y are con- 



sidered known functions of r, then the evolution equation 
for y takes the form 



d d 



(13) 



where the function S(y,r) contains no derivatives of y. 
If we perturb the quantity y about its time-independent 
solution by writing y — > y + £, then jl3| ) yields, to first 
order in £, 



(14) 



where R(r) does not depend on £. 

For each of the ER evolution equations (0) there is 
a corresponding perturbation equation of the form (|l4|) . 
Each perturbation equation has a different function R(r) 
that depends on the right-hand side of the corresponding 
evolution equation. We will see that the form of R(r) is 
what determines whether a particular evolution equation 
is individually stable. 

For the simple case in which (3(r) and R(r) are con- 
stants and /? > 0, the solution to (|l4|) on r G [2M, oo] 



£(r,t)=6>(r + /3i)e 



Rt 



(15) 



where £o( r ) is the initial perturbation at t = 0. The 
stability is determined by the sign of R: If R > (as- 
suming that the initial perturbation falls off with radius 
more slowly than e~ rR ^), the perturbation grows expo- 
nentially with time; if R < (assuming that the initial 
perturbation grows with radius more slowly than e r l fl l/ /3 ), 
the perturbation decays. 

For the more realistic case of nonconstant R and /3, 
the solution to ( O ) is more complicated than ([15]) and is 
considered in Appendix |A|. Nevertheless, one can roughly 
determine whether a given ER evolution equation is in- 
dividually stable by examining the sign of the function 
R{r) associated with that evolution equation. 

Applying this criterion to the ER evolution equa- 
tions (0), we find that R(r) is everywhere negative for all 
but four of these equations, indicating that these equa- 
tions should be stable to small perturbations. The four 
remaining equations have positive R(r), suggesting that 
they might be unstable. If R(r)[ y ] denotes the function 
R{r) associated with perturbations of the variable y, then 
the four positive R(r)[ y ] are 



R(r) [Kt] =ANK t = 
R(r) [aT] = 2NK T = 
R( r )[M rl 



2z 3 



M(l + z)(l + z 2 )' 
z 3 

M(l + z)(l + z 2 )' 
d 



ANK T + —P 
or 

z 3 (2 + 3z + 4z 2 + 5z 3 ) 
2M(l + z) 2 (l + z 2 ) 2 



(16a) 
(16b) 

(16c) 
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R(r) [Lrr] = N(AK T - 5KI) + 2— (3 

z 3 (20+ 19.Z + 18z 2 + 17z 3 



4M(l + z) 2 (l + z 2 ) 2 



(16d) 



where z = 2M/r and the expressions in terms of z have 
been obtained from the analytic solution (p"l|). 

We can test whether perturbations of individual evo- 
lution equations are indeed unstable by modifying our 
code so that a single dynamical variable may be evolved 
in time while all other variables, including the shift, are 
held fixed to the analytic solution ( flj|) . We find numer- 
ically that all evolution equations (0) are individually 
stable except (f7nj), the equation for L rr . 

Our above analysis predicted that the L rr equation 
should be individually unstable because it is associated 
with a positive R(r). However, it also predicted that the 
Kt, ax, and M t t equations should be unstable for the 
same reason. As shown by a more detailed analysis in 
Appendix |A|, the Kt, &t, and M t t equations are sta- 
ble because their corresponding values of R(r) are much 
smaller in magnitude than the value of R(r) associated 
with the L rr equation. 

The growing mode allowed by the L rr evolution equa- 
tion (7n) can be described as a continuum instability: it 
depends only on the equation itself and the equilibrium 
solution, and not on numerics. The only role of numer- 
ics is to produce the initial perturbations that excite the 
unstable mode. 



V. MODIFIED EVOLUTION EQUATIONS 
A. Modifications for stability 

The large positive R(r) associated with perturbations 
of L rr originates from the term N L tt {AKt — 5K^) that 
appears on the right-hand side of the L rr evolution equa- 
tion ( fFn| ) . This term must be modified if the L rr evolu- 
tion equation is to be made individually stable. There 
are several ways to accomplish this. 

One possibility is to change variables. If one evolves 
some quantity QL rr instead of L rr , where Q is some com- 
bination of the other ER variables, then perturbations of 
QL rr will be governed by (|l4| ) with some new value of 
R(r). By careful choice of Q one hopes to obtain a more 
favorable (more negative) R(r). For example, the evo- 
lution equation for the quantity B 2 r 2 L rr yields R(r) = 
N(2K T - 5K~Z) + 2d(3/dr, which is still positive but is 
slightly smaller in magnitude than ( |l 6d| ) . Similarly, the 
evolution equation for U r yields R(r) = N(2Kt — 3if, r ). 
However, there are two reasons why such a procedure 
is unattractive as the sole method of stabilizing the L rr 
equation. First, the ER equations are linear in Lij, M^, 
dij, and aoi (but nonlinear in the other variables), and 
evolving QL rr where Q is anything other than the metric 
functions or the lapse would spoil this linearity. Second, 



in order to make R(r) nonpositive everywhere by evolv- 
ing the quantity B n r n L rr /A m , it turns out that the re- 
quired value of n is large enough that B n r n L rr /A m grows 
with r, hampering our ability to impose an accurate outer 
boundary condition. 

Another approach is to use the constraint equations 
to eliminate the troublesome term that appears on the 
right-hand side of the L rr evolution equation (7n). In 
order to avoid changing the hyperbolic character of the 
evolution equations, one must use only constraint equa- 
tions that are algebraic, that is, those that contain no 
derivatives. Fortunately, many of the ER constraints are 
algebraic. For some constraints this is merely a result of 
spherical symmetry, but several ER constraint equations 
are algebraic even in the general case of three spatial di- 
mensions plus time. In spherical symmetry, the algebraic 
constraints are Eqs. (jBc|), (|8d|), (|8e|) , (S_h), and (|J). An 
additional algebraic constraint can be formed from ( |8a| ) 
and (8b) by eliminating the derivative of T r T, yielding. 

Lrr ~ ~r2 ~ 



2Li 



A 2 



K 



2ai 



A 2 



(«,. 







(17) 



Because we wish to modify the L rr term on the right- 
hand side of (f7n|) for the case in which all variables ex- 
cept L rr are fixed to the analytic solution, the only rel- 
evant algebraic constraints are those that involve L rr , 
namely Q and ([P?]). 

We have found several methods of obtaining an indi- 
vidually stable evolution of L rr - These all involve the 
use of algebraic constraint equations, and some also em- 
ploy a change of variables. We have had the most success 
with the following approach: First eliminate Lt from ( [Sc) ) 
and (|l7]) to obtain 



AA 2 



4B 2 r 2 



A 2 (K 7 



K 



r 2 \ 



o. 



(18) 



Then write down the evolution equation for the quantity 
L r r , and add N(AK T - 5K^)/A 2 times Eq. (|ij|) to the 
right-hand side, yielding 



— —M r 

A 2 dr rr 



N 
3a r — T 



2K r r v; 



TrT 

Br 



Mr 



A 2 



M Tr T rT 
A 2 Br 



-(5K r r - 4K T ) 



1 - 



B 2 r 2 

2 



L r T 
AA 2 



+ ^(6K T + 4K r r ) + -^{2K T + m r r ) 
+K T {2Kf + 3K T K; - AK T 2 ) 



(19a) 



Because we now evolve U r instead of L rr , we also choose 
to evolve M^ r instead of M rrr . This preserves the sym- 
metry between the L-M pairs of evolution equations that 
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make up wave equations. The evolution equation for M r r 
is 



d M r r 



-N 



a r Ll + 2Kfa r + AK r r M, 



r M r 

r rr 



(19b) 



Evolving M£ r has an ad ditio nal advantage: perturba- 
tions of M£ r governed by (19b) have a smaller (more neg- 
ative) R(r) then perturbations of M rrr governed by (|7r. 
so perturbations of M£ r should decay more rapidly. 



B. Results 




FIG. 6. The £2 norm of the error in A versus time, 
computed for three resolutions using the modified evolution 
equations. The outer boundary is at r max = 128M and 
Ai/Ar = 3/4. For t > 5M the growth is only linear in 
time, and the code runs much longer than for the case shown 
in Figure bl. 



Figures |6j — |S| show the £2 norms of the error in A, the 
error in L rr , and the Hamiltonian constraint for simu- 
lations in which we sol ye t he modified evolution equa- 
tions ( |l9| ) in place of (7m) and (^n|). The numerical 
method used in these simulations is identical to the one 
used to integrate the unmodified evolution equations in 
section [II. We use a larger outer boundary radius, 
r max = 128M, to suppress outer boundary difficulties 
that become important at late times. 

For the same grid resolution, our code integrates sev- 
eral orders of magnitude farther in time when using the 
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FIG. 7. The £2 norm of the error in L rr versus time, shown 
for the same three cases as in Figure ^[ 



modified evolution equations than when using the un- 
modified ones. The large errors that grow on a time scale 
of 10 — 100M in Figures are not present in Figures 

Instead, numerical errors increase linearly with time 
(or slower than linearly) for over 10 000A/ until difficul- 
ties associated with our treatment of the outer boundary 
eventually halt the simulation. 

The errors in all dynamical variables except N and 
Br exhibit the same linear growth as seen in Figures ^ 
and [?]. Errors in N and Br are instead dominated by 
outer-boundary effects that grow rapidly and eventually 
terminate our code. Figure shows the error in the lapse 
function N at various times, plotted as a function of ra- 
dius for several simulations with different outer bound- 
ary radii r max but with the same grid resolution Ar and 
time discretization At. Increasing the outer boundary 
radius suppresses the rapid growth of outer-boundary- 
related errors at late times and allows for much longer 
simulations. It should also be possible to improve our 
results by modifying our outer boundary condition, but 
the integration times achieved by our code are already 
beyond what should be necessary for modeling interest- 
ing 3D astrophysical problems such as black-hole binary 
coalescence. 



VI. DISCUSSION 

The success of our free evolution scheme when solv- 
ing the modified ER equations is strong evidence that 
the growing continuum mode identified in section IV is 
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FIG. 8. The £2 norm of the Hamiltonian constraint versus 
time, shown for the same three cases as in Figure ^[ There is 
no significant growth at late times. 



responsible for the instability discussed in section III C . 



The key modification required to suppress the instabil- 
ity was the removal of a term on the right-hand side 
of the L rr equation, the very term that our analysis 
in section [V singled out as problematic. Although we 



have also improved the performance of our code by using 
I7 r and M£ r as dynamical variables instead of L rr and 
M rrr , we have verified that making this change of vari- 
ables alone, without removing the troublesome term in 
the L r r equation by means of constraints, yields r esults 
only marginally better than those shown in section III C . 
Conversely, removing the unstable term and evolving L rr 
and M rrr instead of L r r and M£ r still allows evolutions 
to thousands of M. 

It is no surprise that t he de tailed behavior of the in- 
stability shown in section HI C is much more complicated 



than that predicted by our simple analysis in section IV 
and in Appendix [A| We considered the evolution of a sin- 
gle variable according to a single linear advective equa- 
tion that possesses only fixed, ingoing (for f3 > 0) char- 
acteristic curves. The ER system is actually a coupled 
system of nonlinear advective and wave equations, and 
its three families of characteristic curves (along the in- 
going and outgoing light cone, and along the normal to 
the foliation of time slices) depend on the solution. One 
could do better than our treatment in section IV by lin- 
earizing (0) about the analytic solution and solving the 
entire system of coupled linear partial differential equa- 
tions; however, our approach is far simpler and appears 
to give the correct qualitative results. 




FIG. 9. The absolute value of the error in N as a function 
of radius, shown at various times for several cases of differing 
r max . All plots have Ar/M = 1/32 and At/Ar = 3/4. The 
simulation with r max = 64M crashes at 12 OOOAf. 



We emphasize that the results presented in section V B 



were obtained using free evolution, and that no con- 
straints have been enforced. Furthermore, we note that 
the modifications discussed in section [v] do not alter the 
hyperbolic character of the system. A different version 
of our code evolves (fjj) while enforcing several algebraic 
constraint equations (specifically, we solve (^d|) for M r r, 
( |j~7| ) and ( j3q ) for Lt, and ( |8ej ) for ao r after every time 
step), and yields evolutions accurate for times on the or- 
der of 1000M. While constraint enforcement allows our 
simulations to remain accurate for far longer times than 
with free evolution of the unmodified ER equations (Q), 
our partially-constrained method eventually succumbs to 
an instability slightly after 1000M. The details of exactly 
how constraint enforcement suppresses the continuum in- 
stability found in section ^ are unknown. 

We have concentrated on a case in which the analytic 
solution is manifestly time- independent, namely, when 
initial data given by Eqs. ([ll]) are evolved using a har- 
monic time coordinate. However, modifying our evolu- 
tion equations also dramatically improves our numeri- 
cal results when initial data are chosen on a minimally- 
modified ingoing Eddington-Finkelstein (MMIEF) jg] 
time slice, so that subsequent evolution using harmonic 
time slicing yields a time-dependent result. Using our 
partially-constrained code, we have shown [^8| that the 
evolution of MMIEF initial data using harmonic time 
slicing relaxes to the solution ( [TT| ) at late times. The 
same result holds for free evolution of the modified ER 
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equations 

It is straightforward 
tion 



to extend the analysis in sec- 



[V to the three-dimensional ER system. In this case, 



it is useful to include couplings between tensor compo- 
nents. For example, three-dimensional perturbations of 
Lij, with all other quantities held constant, obey 



d_ 

dt 



(20) 



which is similar to (Uj) except that here £ is a column 



vector containing [8L XX ,8L X 



,5L ZZ ) and R(x,y,z) 



is a matrix. For perturbations about the spherically- 
symmetric solution (pd|), we find that the largest eigen- 
value of R(x, y, z) is given by the same expression (16d) 
as in the spherically symmetric case, indicating that the 
three-dimensional ER equations should suffer from the 
same instability as their spherically-symmetric counter- 
parts. Applying the same analysis to the Kij evolution 
equation in the ADM system (using the same analytic 
solution (pd])) results i n eig envalues of R(x, y, z) that are 
the same size as Eq. (16b) and applying it to the £«, 



Hki.j, K^, and T^- equations in the Einstein-Bianchi sys- 
tem yields eigenvalues of R(x, y, z) that are no larger 
than 3/2 the size of Eq. (16b), so we expect that the type 



of continuum instability we find in the ER system should 
not be present in cither of these two other formalisms. 

Although our stability analysis makes use of the an- 
alytic solution (11), in principle any other solution can 
be used instead as a background for perturbations. Be- 
cause the form of the ER evolution equations given 
by (0) assumes harmonic slicing, the only relevant time- 
independent solution is (0). However, for the case of 
the Einstein-Bianchi or ADM system evolved using a 
different gauge, one might be interested in a different 
background solution. The features of the background so- 
lution that are important for determining stability are 
the signs and relative magnitudes of components of K%j 
and derivatives of (3 Z . We note that these features are 
approximately the same for the Schwarzschild solution 
on time-independent MMIEF slices as they are for the 
Schwarzschild solution on time-independent harmonic 
slices, so one obtains similar stability criteria in both 
cases. 

In the case of the ER equations, we are fortunate to 
have algebraic constraints that can be used to modify 
the evolution equations without affecting the hyperbolic 
character of the system, even in three dimensions. How- 
ever, not all the ER constraints are algebraic, and it is 
unclear in the three-dimensional case which constraints 
must be used in order to suppress instabilities. In par- 
ticular, Eq. dl8|), which seems necessary for removing 
the growing mode, is not algebraic in three dimensions. 
This is because Eq. ( p"8| ) results from eliminating second 
derivatives of the metric from Eqs. ( |Sa| ) and (pb]); the 
three-dimensional equivalent is forming a linear combina- 
tion of components of Eq. (Q) that eliminates all second 
derivatives of gij appearing in the Ricci tensor R\j , and 
is not possible for a general spacetime. 



One might ask why we do not use ( |8c[ ) instead of (n| ) 
to obtain a stable evolution equation for L r r , since (|k) 
is algebraic in the general three-dimensional case. The 
answer is that it is possible to use (|8^) to obtain an indi- 
vidually stable evolution equation for L r r . However, doing 
so introduces a term containing Lt on the right-hand side 
of the L r r evolution equation, where no such term existed 
previously. This term generates a continuum instability 
in the coupled L^-Lt system (where all variables except 
U r and Lt arc held fixed to the analytic solution). 

To better understand why (|8c]) alone cannot stabilize 
the ER equations, consider as fundamental variables not 
L rr and Lt, but instead the trace and the trace- free 
part of L^, which in spherical symmetry are given by 
L\ = L r r + 2Lt and L TF = L r r — Lt- If one constructs 
evolution equations for L\ and L TF , one finds that per- 
turbations of L TF , holding L\ and all other ER variables 
fixed, obey (Oh with 



R(r) 



N(10K T - 7K r r ) _ z 3 (48 + 41 z + 34z 2 + 27z 3 ) 



12M(1 



! (1 



(21) 



The perturbations grow rapidly with time because R(r) 
is large and positive. The source of the problem is a 
large, positive L TF term on the right-hand side of the L TF 
evolution equation. Because ( pc| ) involves only the trace 
of L^ and not its trace-free part, this equation cannot be 
used to eliminate the L TF term and thus cannot be used 
to stabilize the system. 

If one wishes to use the ER formulation in a 3D free 
evolution, one must find a way of dealing with the un- 
stable continuum mode afflicting the ER evolution equa- 
tions. Unfortunately, the above analysis suggests that in 
3D, this cannot be done in a simple way using algebraic 
constraint equations. Accordingly, for 3D simulations it 
may be more fruitful to pursue other hyperbolic formu- 
lations such as the Einstein-Bianchi system, which, ac- 
cording to our analysis, should not suffer from this type 
of instability. 
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APPENDIX A: 



1. Solution of Eq. ( |14| ) on infinite domain 

Solutions to (jlj) propagate along characteristic curves 
r = r(t) that depend only on the shift vector and are 
defined by 



(Al) 



Each spacetime point (r, t) intersects exactly one charac- 
teristic curve. If we define s(r, t) to be the radial coor- 
dinate at which the characteristic curve passing through 
(r, t ) int ersects the initial slice t = 0, then for /3(r) given 
by ( lid ) we can integrate (|Al[ ) to find a relation between 
s, t, and r: 



t 

2M 



ln£ 
r 



1 / r \ 3 
3 V2M7 

if—) 2 

2 \2M/ 



(;)' 
(;)" 



r 
2M 



s 

- - 1 

r 



(A2) 



Treating s and i as independent variables, we can write 
Eq. ( |l4| ) in the form 



dt 



(A3) 



Each value of R{r) listed in (|16| ) can be written in the 
form 



R(r) 



z 3 (a + bz + cz 2 + (b + c — ajz 3 ) 
2M(l + z) 2 (l + z 2 ) 2 : 



(A4) 



where a, b, and c are constants and z = 2M/ r. Using 
this expression for R(r) we can integrate ( |A3| ) together 
with ( |A2|) to obtain the general solution 



£(r,t) =£„(*) - 



1 + 2M/r 
1 + 2Af/s 

' 1 + (2A//r) 2 
1 + (2M/s) 2 



(c-o)/2 



(A5) 



where £o( r ) denotes £ on the initial slice t = 0. 

For a fixed value of r we have s r at late times, 



so (A2) reduces to t ~ s 3 /12M 2 and (A5) reduces to 



£(r,t)~£o(12 1/3 A/ 2/3 t 1/3 ) 



12Af 2 t 



a/3 



(A6) 



where time-independent factors have been dropped. If £ 
initially falls off like r~ m , then for a fixed r it will behave 
like t( a_m )/ 3 at late times. For a > m + 3, perturbations 
will grow superlinearly with time, but for a < m + 3 the 



growth is at most linear (for a < m the perturbation is 
actually damped), so it does not represent an instability. 
For the L rr equation (a = 10) to be individually stable, 



numerical errors must fall off at least as fast as 



For 



the Kt equation (a = 4) to be stable, the leading-order 



falloff rate must be no slower than 



The AI r T and ay 



equations (a = 2) will be stable even if numerical errors 
grow with radius, as long as these errors grow no faster 
than r. 

Empirically, we find that the dominant numerical er- 
rors in the wavelike variables (L rr , Lt, M r T, a rr , and 
ciQ r ) fall off like r _1 and propagate outward from the 
strong-field region near the hole. This is what one would 
expect for modes that behave like gravitational radiation 
(these modes are not allowed in spherical symmetry but 
nevertheless can be present in numerical error terms). 
The dominant errors in other variables also propagate 
outward from the strong- field region, and fall off either 
like r" 1 or r~ 2 . These falloff rates explain our observa- 
tion that the L rr equation is individually unstable but 
the Kt, M t t, and ax equations are individually stable. 

For background solutions other than (11), the forms 



of /3(r) and R{r) will be different, so the details of the 
solution ( |A5| ) will change. For example, if one takes the 
MMIEF solution as a background (this is not relevant for 
Eqs. (0) because the MMIEF solution is not preserved un- 
der harmonic slicing, but is relevant for other systems of 
evolution equations to which one might apply this anal- 
ysis), R(r) typically falls off like r~ 2 instead of r -3 , and 
/3(r) falls off like r _1 instead of r~ 2 , so the stability crite- 
rion becomes a < m+2 instead of a < m+3. At the same 
time, the coefficient a is typically smaller for the MMIEF 
background, so both the MMIEF background and the 
background ([ll]) yield similar predictions for stability. 

Furthermore, note that our stability criterion can be 
applied to coupled evolution equations as long as there 
are no couplings through derivatives. For example, con- 
sider the coupled system consisting of all ER variables ex- 
cept M rrr , M r T, and ao r . If M rrr , M r T, and a^ r are held 
fixed, the perturbation equation for the thirteen other 
variables can be written in the form (]lj), where in this 
case £ is a thirteen-element column vector and R(r) is 
a 13 x 13 matrix. To determine stability, one examines 
each eigenmode of the perturbation equation in the man- 
ner described above. An example in which this analysis 
cannot be used without modification is the coupled sys- 
tem consisting of L rr and M rrr . In this cas e, the spatial 
derivatives of L rr in the M rrr equation (7m) and the spa- 
tial derivatives of M rrr in the L rr equation (7n) prevent 
one from writing down a matrix perturbation equation of 
the form (|l4|). Instead, the perturbation equations pos- 
sess more than one family of characteristic curves, so the 
solution is more complicated. 
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2. Solution of Eq. ( |l4| ) on finite domain 

In numerical simulations one often does not have a 
domain that extends to r — oo, but instead one imposes 
an artificial boundary condition at some finite value of r. 
For simplicity, consider a Dirichlet condition: assume £ 
is fixed to some constant value £o at the outer boundary 
r = rQ. If we let io( r ) be the time it takes for information 
to propagate from the outer boundary ro to some radius 
r < ro, then for (i, r) such that t > to(r), the solution 
of (0) is time-independent, and for (i, r) such that t < 
to(r), the sol utio n is the same as for the case considered 



in Appendix A 1 . 

For /3(r) given by (lid), the time it takes for informa- 
tion to propagate inward from radius s to radius r < s is 
given by (|A2|). In this case, for R(r) given by (A4) the 
time-independent solution is 



1 



- 2M/r \ 
f 2M/r ) 

1 + (2M/r) 2 
1 + (2M/r„) : 



(A7) 



For ro 3> r, the time-independent solution behaves 
roughly like r~ a . 

One consequence of the above analysis is that if one 
uses a Dirichlet outer boundary condition and an unsta- 
ble mode of this type is present (that is, if numerical 
perturbations fall off more slowly than r~ a ), then the in- 
stability will become more severe if the outer boundary 
location is moved to a larger radius. This is because the 
unstable mode has more time to grow before the time- 
independent state is reached. We have verified this nu- 
merically for the simple case of the L rr evolution equation 
solved with all other variables held constant. 
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